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Random time averaged diffusivities for Levy walks 
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Abstract. We investigate a Levy- Walk alternating between velocities ±iio with opposite sign. The sojourn 
time probability distribution at large times is a power law lacking its mean or second moment. The first case 
corresponds to a ballistic regime where the ensemble averaged mean squared displacement (MSD) at large 
times is <^a;^) ex i^, the latter to enhanced diffusion with (s^) oc f", 1 < i^ < 2. The correlation function 
and the time averaged MSD are calculated. In the ballistic case, the deviations of the time averaged MSD 
from a purely ballistic behavior are shown to be distributed according to a Mittag-Leffier density function. 
In the enhanced diffusion regime, the fluctuations of the time averages MSD vanish at large times, yet 
very slowly. In both cases we quantify the discrepancy between the time averaged and ensemble averaged 
MSDs. 

PACS. 5.40.Fb, 02.50.-r 



1 Introduction 

Dispersion of Brownian particles is described on the macro- 
scopic level by Pick's second law which states that the 
local change in particle density is proportional to the neg- 
ative gradient of the local particle flux, where the diffu- 
sion constant D is the proportionality factor. Accordingly, 
the density of the particles is a spreading Gaussian with 
the mean squared displacement (MSD) going linearly with 
time, (x^) = 2Dt. Einstein Qj established a relationship 
between the macroscopic quantity D = {{dxY) /2{t) and 
the underlying stochastic process with jump lengths dx 
of variance {{dx)"^) and the average time passing between 
two jumps (r). For a single Brownian particle, the time 
averaged MSD (TAMSD) P has to be considered. 



P = 



1 



t-A 



t-A 



[x{t' + A) - x{t')f dt'. 



(1) 



Here the averaging was performed over all displacements 
along the single trajectory occuring during a fixed lag time 
A within the measurement time t. Due to the station- 
ary increments, for a Brownian motion the time averaged 
MSD will attain the limit 5'^ = 2D A for large times t. 
Brownian motion is therefore ergodic, ensemble and time 
averages are equal, (x^) — S'^. It is important to note that 
in practice an estimation of ensemble averaged diffusion 
constants from single particle trajectories even for normal 
diffusion is a subtle issue due to the lack of statistics. How- 
ever one can exploit the ergodic property of such processes 
in order to find suitable estimators [I], [S]- 

Unlike in normal diffusion, in strongly disordered me- 
dia the mean squared displacement often does not grow 



linearly with time, but according to a power law (a;^) oc f^ . 
The case < i^ < I corresponds to subdiffusion, I < i/ < 2 
indicates enhanced diffusion or superdiffusion. In disor- 
dered and complex systems time averages can differ con- 
siderably from the corresponding ensemble averages. The 
increasing employment of single particle tracking tech- 
niques makes it indispensable to understand these differ- 
ences in order to interpret the experiments and gain in- 
sight to the mechanisms underlying anomalous transport. 



Examples for enhanced diffusion are experiments on 
active transport of microspheres 4 , of polymeric particles 
[5] or of pigment organelles |6J in living cells. Enhanced 
diffusion in living cells is promoted by molecular motors 
that move along microtubules or the cytoskeleton [7], [5]. 
In vivo experiments usually examine the trajectories x{t) 
of single particles and hence assess the TAMSD Eq. ([T]) in- 
stead of ensemble averages. Measurements often find this 
quantity to be a random variable, so that ensemble average 
and (single trajectory) time average differ. For this behav- 
ior several reasons come into question, either variations in 
the probed environments or cells, ergodicity breaking or 
too short measurement times. In Refs. [1], [5] the observed 
enhanced diffusion was described within the framework of 
generalized Langevin equations (GLE), which is Gaussian 
and ergodic [5]. Using this theory, the experimentally ob- 
served exponents characterizing the anomalous transport 
were reproduced. However, this Gaussian approach can- 
not explain the multiscaling of moments found for the en- 
hanced motion of polymeric particles in living cells [5] , a 
feature that seems to be more consistent with a Levy walk 
scheme. 
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Levy flights describe enhanced diffusion in terms of 
random walk processes where the distribution of the par- 
ticle displacements lack the second or even the first mo- 
ment (i.e. give rise to a Levy statistics) . Levy flights have 
been used in the past to describe phenomena as diverse 
as the dispersal of bank notes [TU] , tracer diffusion in sys- 
tems of breakable elongated micelles [H] or animal forag- 
ing patterns [I^- However, Levy Flight models can lead to 
unphysical behavior with regard to velocities since Levy 
flights are characterized by extremely large jump lengths, 
corresponding to the heavy tails in the jump length distri- 
butions. Levy walks address finite velocities by either pe- 
nalizing long instantaneous jumps with long resting times 
(jump models), or by letting the particles move at a cer- 
tain velocity for a certain time or displacement, and choos- 
ing a new direction (or velocity) and sojourn time accord- 
ing to given probabilities (velocity models) [13], [13]. 

The theory of Levy walks finds a wide range of ap- 
plications. Experiments with passive tracer particles in a 
laminar flow have shown that the flight times and hence 
displacements within the resultant chaotic trajectories of 
the tracer particles can exhibit power-law distributions 
|15| . Likewise, the motion of tracers in turbulent flows 
can be described as Levy walks [TB]. Another example 
is the stochastic description of on- off-times in blinking 
quantum dots, where the intensity corresponds to the ve- 
locity of a particle alternating between states of zero and 
constant nonzero velocity. Nonergodic behavior was found 
in the correlation functions for sojourn time distributions 
lacking their mean jI7],[T2]. Another example is the dy- 
namics of cold atoms in optical traps [Hj and the re- 
lated Brownian motion in shallow logarithmic potentials 
[20] . or perturbation spreading in many-particle systems 
[51]. Moreover, also deterministic systems such as certain 
classes of iterated nonlinear maps may show enhanced dif- 
fusion. The chaotic behavior of resistively shunted Joseph- 
son junctions manifests itself in an anomalous (determin- 
istic) phase diffusion, which can therefore be modelled by 
means of such maps |22j . In turn, the enhanced diffusion 
behavior emerging from such iterated maps can be mod- 
elled stochastically using the Levy walk approach (in par- 
ticular velocity models) [15] . 

In this article we study Levy walks in one dimension 
where the persistence times in the positive- or negative ve- 
locity state are drawn according to a probability density 
function iP{t) with either first or second moment lacking. 
The first case is referred to as the ballistic case, the latter 
as the subballistic or enhanced case. Recently, the dynam- 
ics emerging from nonlinear map similar to the Levy walk 
was investigated [22] • However, this work did not address 
the fluctuations of the TAMSD. Fluctuations were consid- 
ered in a very recent numerical study for the special case 
of Levy walks exhibiting enhanced diffusion pT , and nu- 
merically and analytically for enhanced and ballistic case 
in a brief publication of the authors [3S] . 

The article is divided into four parts. The first one is 
dedicated to the ballistic case of a Levy walk, the second 
one to a Levy flight with a step size distribution lacking 
the second moment, and the third one to the enhanced 



Levy walk case. For both the ballistic and the enhanced 
case we first review briefiy occupation times, propagators 
and ensemble averaged mean squared displacements. Then 
we turn to the ensemble averaged quantities such as cor- 
relation functions and ensemble averages of the TAMSDs. 
Finally, we investigate the distributions and properties of 
the fiuctuations of the TAMSDs. In the second part, we 
investigate for comparison the Levy fiight corresponding 
to the enhanced case. We also provide the fluctuations 
of time averages of Levy flights, using simple arguments, 
thus adding to the work in 35 who addressed this prob- 
lem rigorously and more generally. 



2 Ballistic regime 

We consider a one-dimensional motion of a particle with 
a two-state-velocity ±vq where the sojourn times in the 
states are drawn from a probability density function (PDF) 
iP{t). Hence the particle has a velocity +vo for period ri 
drawn from iP{t), after that switches to velocity —vq and 
remains in this state for another period T2 also drawn from 
'0(r). This process is then renewed. In particular, this PDF 
is chosen such that it lacks its first moment with a power- 
law decay at large times, iP{t) ^ A/r{—a)T~'^~°' with 
< a < 1. Particularly in the simulations we will use 



Hr) 



-l-a 



T> 1 

else . 



(2) 



The Laplace transform of Eq. ^ in the small-w-limit is 



^{u) ~ 1 - Au" 



(3) 



with A — r(l — a) and u being the Laplace conjugate of 
r. This relation can easily be derived via Laplace transfor- 
mation of L \p{T')dT' ~ 1 — T^", which yields ^V'('") — 
i (1 — r{\ — ck)u") by using convolution and Tauberian 
theorems. 

In the limit of long times t the particle moves ballisti- 
cally so that the mean-squared displacement is (x^) = (1 — 
a)t^ [33] , [37] . Recently, a similar system had been gener- 
ated in the context of deterministic superdiffusion and the 
ensemble average of the time averaged mean squared dis- 
placement ((5^) was derived [33]. Processes whose temporal 
dynamics is governed by heavy-tailed PDFs lacking their 
mean exhibit ageing [28] which manifests itself in a de- 
viation of the ensemble averaged mean-squared displace- 
ment from the TAMSDs which are random themselves. 
Therefore, the fiuctuations of TAMSDs are an important 
signature of this process. 



2.1 Occupation fraction and particle position 

The distribution for the fraction of time z± = t±/t spent 
in one state (positive or negative velocity) after a large 
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time t is given by Lamperti's law |29j . 



Pocciz±) 



smTTQ; 



-^(i--±r 



^2q 



(1 — z±) + 2c0S7rQ;zJ (1 



z±r 

(4) 
a generalization of the arcsine law which is reproduced for 
the case a ~ 1/2. The particle position x{t) is given by the 

integral over the velocities /„ v{t')dt', and the temporal 
mean of the velocity is x/t = fo(i+ — t-)/t- The distribu- 
tion of this quantity in the limit of large times had been 
calculated in the context of the mean magnetization of a 
two-state system |30] . A related problem is the calculation 
of the integral over the intensity of emitted light in blink- 
ing quantum dots. In that case one state corresponds to 
the "on" -state, i.e. the emitting state, and the other state 
to the "off" -state where no light is emitted [TT. 

The probability to find the particle at position x at 
time t for large times can be obtained in terms of the 
scaling variable z = x/{vot) by using Eq. (|3]) and change 



of variables. We have 

dz 



2z, 



1, hence 



Piz) 



p{z) = 



Po 



(z+), which finally results in 



i and 



2sin7ra (l — z^) 



TT ( (1 + zf" + (1 - zf" + 2 cos Tra (1 - z^Y 



(5) 

Fig. [T] shows this distribution of the scaled particle posi- 
tion for two different values of a. 




wHm 



-1.0 -0.5 



-1.0 -0.5 



Fig. 1: Distribution of the position of particles at t — 10^ 
for a = 0.5 (left) and a = 0.7 (right panel) in terms of the 
scaling variable z [vq = 1). Sample size 10^, '0(t) is given 
byEq. ©. 



In order to find ((5^), we first derive the Levy walk cor- 
relation function {x{ti)x(t2)), which turns out to exhibit 
ageing. The position correlation function {x{ti)x{t2)) is 
related to the velocity correlation function {v{si)v{s2)) 
via 



fti rt'2 

{x{ti)x{t2)) ^ { I v{si)dsi / v{s2)ds2 
Jo Jo 

dsi / {v{si)v{s2)) ds2 

/O Jsi 

ds2 / {v{si)v{s2)) dsi 







(8) 







where we took into account that ^2 > ^i- Using the ap- 
proach of Ref. [3D] we obtain for the velocity correlation 
function {v{ti)v{t2)) 



){tl)v{t2)) ^vlY,{~^TPn{h,t2) 



(9) 



where Pn(ti,t2) is the probability of the velocity to switch 
its sign n times within the time interval [^1,^2], i-e. for 
even n between ii and ^2 we have v{ti)v{t2) = Wq, and for 
odd n, v{ti)v{t2) = —Vq. 

In the scaling limit where ti and ^2 are large and Eq. 
(|31) applies, the particle gets stuck in the -I- or — state 
for times of the order of the measurement time due to the 
lacking first moment of the sojourn time distribution tpir) . 
Therefore, only the first term n = is relevant in Eq. ([^. 
The corresponding probability ^0(^1:^2) of the velocity not 
switching its sign from a given ii on up to ^2 is called the 
persistence probability, to which the velocity correlation 
function is proportional, {v{ti)v{t2)) — VQpo{ti,t2). Let 
us denote the first waiting time from an arbitrary time 
ti up to the next switching event by r/, see Fig. [2] This 
first waiting time is called the forward recurrence time, 
and its PDF -0j ^^(t/) differs from ^(r) since ii does not 
necessarily coincide with a renewal event. 



2.2 Ensemble average of S^ 

First we will analyze the ensemble averaged TAMSD: 
1 



rt-A 



{S') = 



t-A 



[x{t' + A) - x{t')Y dt' 



U ^f 



T/'+l T/+2 



iQ\ Fig. 2: Sketch of the forward recurrence time Tf. Ticks 
symbolize renewal events. 



Changing the order of integration and ensemble averaging 

in Eq. ([S]), we get With ^2 —ti+ A, we express the persistence probabil- 

* A n o ity as (30| 

^'-^ {x'if + A)) + {xHt')) 2{x{t'){t' + A)) ^ , '-' 



(<52) = 



t-A 



-dt' 

(7) 



Poih,h+A)= i;fj,{Tf)dTf. (10) 

Ja 
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In terms of the scaling variable 9 = A/ti, the pdf of the recalling that ti/t2 < 1 and taking the limit a -> we are 

forward recurrence time to take a value of A at time ti left with 

reads in the scaling limit where also A is large: „ 

lirria^o {x{ti)x{t2)) = Wotii2 

r ; io\ — sinTra 1 

lim "tpfyu] — . {II) where we used de I'Hospital's rule. Simulations of the sys- 

tem for moderate a show a good agreement with theory 
The limit theorem for forward recurrence times V'/.ti Eq. Eq. ([T3|), see Fig. [31 
(fTTj) is due to Dynkin [31]. Hence, 



{v{ti)v{h + A)) =vl 



siuTTO: 


1 


d9 


IT 9^ 


'{1 + 9) 


sinTTO; 


1 





^1^.(1 _^)a 



dC 



, sm IT a 



S(^;a,l 



(12) 



where B{y;a,b) — J^ duu"-^^{l — u)''^^ denotes the in- 
complete Beta- function [32] . Note that this expression yields 
only real values for ^2 > ii- In the case of ii > t2, the 
ii and ^2 in Eq. ([T^ have to be interchanged. Inserting 
Ea. (fT2)) into Eq.® and using integration by parts we find 

h 




1 x10" 



{Xih)x{t2)) = V^ 



-tiB 



tit2B — ;a, 1 



t2 



Fig. 3: Simulational results for {x{ti)x{t2)) for a = 0.5 
(upper) and a — 0.7 (lower graph), and the respective 
theoretical predictions (solid lines), Eq. p^ . i2 was fixed 
at 10*, vo = 1. Sample size 10^. 



t2 

h 

t2 



il + a, 1 



:t?B f — ;-l + a, 1-a 



.2 



The asymptotic behavior of the position-autocorrelation 
function Eq. (TT^ is 



-a^t 



(13) {x{h)x{h+A)) 



In particular, for ti = t2 we find the MSD 

{x'{h))^{l-a)vltl, 



7ra(l-Q2) 1 ^ ti ) (^r,\ 



(14) 



in agreement with [55] , [37] , [50] • The theoretical autocor- 
relation functions increase in this case with increasing time 
difference. For normal diffusion we have {x{ti)x{t2)) = 
2Dmh\{tx,t2), so that {x{ti)x{t2)) / {x^ {ti)) = 1 for ^2 > 
ti. In contrast, the behavior of the Levy walk is governed 
by long periods of ballistic motion. Thus, it exhibits strong 
correlations compared to normal diffusion which are due 
to the long sticking times in the positive or negative veloc- 
ity states. In the limiting case a — >■ the particle remains 
in state +wo or — wq throughout the measurement so that 
x{t) = ±VQt with probability 1/2 for either sign. There- 
fore we expect the purely ballistic, deterministic behavior 
{x{ti)x{t2)) — VQtit2 for the position-position correlation 



Note that we made again the transformation from (ti,i2) 
to ti, A = t2 — ti and ti < t2- Inserting the above results 
for the correlation function Eq. (I13p and mean squared 
displacement Eq. P^ into Eq. ([7]), integrating by parts 
and using again the integral definition of the incomplete 
Beta function, we obtain the ensemble averaged TAMSD. 
In the limit A/t <C 1 we get 



m 



A 



2 sin^« 2A^^y 



(16) 



Note that in fact the short time behavior of the correla- 
tion function is not negligible in the integral Eq. ([7]) since 
it affects the long-time behavior of the ensemble-averaged 
TAMSD (S^). It is important to point out that even the 



function, and hence {x{ti)x{t2))/{x^{ti)) = izAi- To see first term of {S^{A)), Eq. dH]), differs from {x^{t)) = 

, , r ^ T^ / f, 1 \ ,• , ^'n(l — ce)t'^ by a factor. This term was also found recently 

this, note that ioi a -^ 0, B { -r;a,l — a] diverges and . -..rr . j.irj. --j.- u Pioi 

' ' V *2 ' y ° ma dirterent context of deterministic maps by [23J . 

the first term in Eq. (|13p is the only term that remains: 



VQtit2 sinvra „ f t 



B 

VQtit2Sm.Tra f ti 



-;a, 1 



Q OO 



r{a + i) 



^-^ r{a) (a + i)i\ \t2 



2.3 Fluctuations of the time averages 

More important are the fluctuations of the TAMSDs, since 
these allow conclusions to be drawn with respect to the er- 
godic properties of the system. Our simulations revealed 
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100 200 



500 



1000 2000 
A 



5000 ixlO'* 



Fig. 4: Deviations from ballistic motion of TAMSDs 
VqA^ — P- versus A for ten different trajectories; a = 0.5, 
vq — 1, t = 10^. Points belonging to the same trajectory 
are connected by a straight line. 



that these fluctuations are quite small compared to the 
value of the TAMSDs, and become even smaller with time 
relative to the ballistic contribution (vqA)'^. The fluctu- 
ations are more pronounced if one looks at the shifted 
TAMSD VqA'^ — (5^, which is the natural random variable 
of this process as will turn out soon. In Fig. |4] we plot 
VqA'^ — S^ versus the lag time A for ten different trajec- 
tories. VqA"^ 



S^ remains visibly random. 



n sojourn time PDFs with the probability of no event after 
the nth one [JT], and is well investigated. Taking i)(u) ~ 
1— Au" in Laplace domain and using the convolution theo- 
rem results inp„(u) = Am"~^ exp [nln(l — Am")]. Laplace 
inversion yields 



Pn(i) 



t 



1 



Q,^l/a^l + lQ 



t 



^l/Q„l/a 



(19) 



la,iit) denotes the one-sided Levy density whose Laplace 
transform is given by exp[— w"] j31] . Hence, with (rit) '^ 
t"{Ar{l + a)) the ensemble average of Eq. P^ becomes 



iS') - vlA^ 



H 



iAr{l + a) 



V^i 



(20) 



which clearly differs from Eq. p^ . Eq. p^ and hence (1^ 
describes a special case of the sojourn time distribution 
Eqs. 0, © fulfilling the relations i^"" > A/{Ar{\ + a)) 
and Zi < 1 so that the first ballistic term in Eq. ([TB]) 
remains larger than second term. In contrast, Eq. (|16p 
requires Z\ ^ 1. From Eq. p5)) we find that the quantity 
{vqA"^ — 5'^) and the amount of switching events within t 
are proportional. 



— 2 A^ 



Therefore, in terms of a new variable 



(21) 



Small A limit: In order to specify the distribution of 
VqA'^ —5-^, we construct the following special case. Con- 
sider again the sojourn time PDF Eq. ([2|) with < a < 1. 
Moreover, let us for now only adhere to small A < 1, so 
that at maximum one change of direction takes place dur- 
ing the interval A. In this case, the TAMSD Eq. ^ can 
be explicitly calculated. Clearly, if there were no changes 
of direction, the TAMSD would be S"^ = v^A^. For a single 
switching event at time ts within a time interval (i, t-\- A), 
the integrand in Eq. ^ becomes 



[x{t + A)-x{t)\ 

r vlA^ for 

\ {2vots - vqA ~ 



t<(ts 

2vQt)^ 



- A) and 
for {ts ~ 



ts <t 

A) <t<t, 



(17) 



Hence, using Eq. ((U, one change of direction within the 
observation time t reduces the TAMSD to a value of S^ = 



v^A^ 



I ^-^'^ for t » Z\ (i — Z\ « t), as shows integration 



of Eq. ^TT} . Two changes of direction result in 5'^ = v'^A'^ — 

2 

2 • I Y"^'^ and so forth. Altogether we find for an amount 
rit of switching events within the observation time t 



t 



Ah--A'nt 
3 



(18) 



where the amount of direction changes rit within (0, t) is 
a random variable. The probability p„(i) of the number of 
events n within (0,t) is determined by the convolution of 



e 



v^A^ - S^ 
\A^ - (P 



nt 
{nt) 



t\ (22) 



and using Eq. ([TO| the rescaled distribution of the TAMSD 
becomes 






r^/"{i + a) 



^l/o 



(23) 



This PDF is the density of the Mittag-Leffler distribu- 
tion, a distribution already encountered in the context of 
TAMSD fluctuations in the subdiffusive continuous time 
random walk [33], [31]. Fig. [S] shows the PDF ^ and 
the respective results for simulations of the Levy- Walk for 
two different values of Q;,_but for large A. It is interest- 
ing to note that {v'^A'^ — S^) differs for small and large A 
regimes, however the distribution of the rescaled variable 
(|22]l does not depend on A. 



Crossover to the large A regirne: The ensemble average 
of the fluctuations of {vqA'^ — S'^) is given by Eq. p^ for 
large 1 ^ Z\ <C i. In the present case where 'tpir) = at 
short times t < 1, the behavior of the ensemble averaged 
TAMSD at Z\ < 1 is given by Eq. i^U^i. This behavior at 
small A constitutes the lower bound for more general iP{t) 
with arbitrary shape at small t. However, the fluctuations 
of the time averages Eq. (|^3|) are governed only by the tail 
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Fig. 5: PDF of C = {{voAf-5^)/{{vQAf-{6^)) at t = 10^ 
for a = 0.5 (left) and a = 0.7 (right panel); A = lO^. The 
histogram shows the result of the simulations, with '0(r) 
Eq. (12), the sohd blue curve the theory Eq. (f^ . Sample 
size 10^, vq = 1. 



times. For this purpose, we calculate the ergodicity break- 
ing (EB) parameter [33] for the shifted TAMSD ^ 



EB 



lim 



{e)~{C}^ _'2r\l + a) 



{£,Y 



r(l + 2a) 



(25) 



where we used Eq. ([25]) . Numerics for a — 0.5 show that 
the EB-parameter for ^ tends indeed to the predicted finite 
value EB = 0.571 (Fig. [7]), i.e. the variable S, remains 
distributed according to Ea. (P^ . 



of the persistence PDF ^(r) and are therefore the same 
for small and large A. Hence we can write 



]A^-6^^X^nu 



x' = 



2v: 



2vf,A-' 

■ 3* 3 
sin-KaA^ 



■na(6-lla+Qa^ -a^)Ar(l+a)t 



Z\<1 
Zi>l 



(24) 



Here x^ gives the deterministic part that governs the en- 
semble mean of the shifted TAMSD, while the full fluctu- 
ations enter via n^, compare Eqs. (1221) . ((23|) . In Fig. [6] we 
plot t{vQA^ — i^"^)) versus A. Simulational results match 
the theoretical short time as well as long time behaviors, 
Eqs. ([20]) and (ITCl) . respectively. The crossover takes place 
in the region of the cutoff of the sojourn time PDF iJj{t) 
at small times, i.e. at Arr ~ 1. 
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Fig. 7: EB-parameter of ^ for a = 0.5. Red circles indicate 
A — 0.5, blue diamonds A = 100, green triangles A = 
5000. The grey solid line indicates the theoretical value 
EB = 0.571, Eq. (125]). Sample size 5000. 
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Fig. 6: t{vlA^ - ((52)) versus A, a ^ 0.5. The small- 
A region is sensitive to the shape of V' (''')■ Lines indicate 
theory Eq. ([^0]) (dashed) matching numerical data (dots) 
at small A, and Eq. ([TO]) (solid) for large A. Note that 
the crossover takes place in the region of the small-time 
cutoff of the sojourn time PDF Eq.(I2]), Acr ~ 1. Sample 
size 10^, vo = 1, t = 10^. 



Finally we demonstrate numerically that the above dis- 
tributions are indeed the limiting distributions at large 



Note that, however, the EB-parameter for the original 
(not shifted) TAMSD S"^ slowly tends to zero for nonzero 
A as ^-2(1-")^ Hence, for the ballistic Levy walk, non- 
ergodicity in the sense of the distribution of time averages 
does not find its expression in the (decaying) fluctuations 
of the TAMSDs 5^ themselves, but in the (persisting) fluc- 
tuations of the shifted and rescaled variable ^. 



3 Levy flight TAMSD 

Before we turn to the behavior of the TAMSD of the Levy 
walk in the enhanced diffusion regime, let us illustrate the 
situation in the related Levy flight. We present a rather 
illustrative than rigorous argument, to avoid complicated 
math. A more general and rigorous treatment can be found 
in Ref. [3S|. The Levy flight is a random walk process 
where at each renewal the displacement x of the walker is 
drawn according to a jump PDF A(a;), but in contrast to 
the normal random walk the jump PDF lacks the second 
moment. With a unit time span passing between consecu- 
tive renewal events, the number of jumps acts as the (dis- 
crete) time variable t. Hence consider the coordinate of a 
Levy flight after t steps as a sum of independent identi- 
cally distributed (i.i.d.) random variables or displacements 
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(y 



asymptotics can be obtained in Laplace domain, 



using the Tauberian theorem: 



Let the {xi\ be distributed according to a two-sided sym- 
metric distribution A(x) = A(— .t) faUing off as a power 
law for large |a;|. 



p{uy) ~ 1 



A 



— ^(1 - o X 

a z 



(32) 



Hence, the sum over these 



A(a:)cx^|:r|-i-" 



1 < a < 2, 



(26) 



= yk in Eq. (1501) is a sum 
over positive i.i.d. random variables distributed according 
to a PDF with a power law tail of exponent —1 — a/2. 
Hence, due to the generalized central limit theorem we 
find in the large t limit [15], [SI]the PDF of C = tS^ j A: 



In particular, for our simulations we used Eq. (1261) with 
Aq = a for \x\ > 1, and X{x) = for \x\ < 1. Then, 
the distribution of the sum Xt will yield a two-sided Levy 
law for large t, La^o{Xt/{tAQr{l — a)/ay^°'), according to 
the generalized central limit theorem [251[?T| . The function 
Lo,aix) is defined as the inverse Fourier transform of 



piO 



1 



{Kk^ty 



-L^.i 



C 



{K^ny 



(33) 



L,^o(fc) = exp[-|fcr] 



(27) 



with k being the Fourier variable [3T1 . This Green function 
of the Levy flight has a similar behavior as the central part 
of the Levy walk Green function when 1 < a < 2, as will 
become obvious in the next section. 

The time-averaged mean squared displacement (TAMSD) 
is defined as [35l 



where K^^ = ^^(1 - f ) and ia/2,i(-) is the one-sided 
Levy PDF given by exp [— u"/^] in Laplace domain [31) . 
A similar result was obtained in [3S] though a slightly 
different scaling was reported. Fig. [8] shows the PDF Eq. 
(P5)) obtained with Mathematica, and the corresponding 
simulational results which perfectly match the theory. 



(52 



1 '-"^ 



1 





(28) 



t-A 
where the integer A is the lag time. We have 

t-A /i+A i+Ai+A 

i—l \ k—i k—i j^k 

The mixed terms X]fc=i X^i^fe ^kXj on average cancel out 
for large enough A, hence we omit them. Moreover we 
assume 1 <?; ^ ^ t so that 



(29) 



t-A 



t-Ai+A 



^— 1 k—i 
t i+A 



t 

d A 



1 ■- ^^^ 

tEE 



^k 



z— 1 fc— i 
fc=l 



We find for the distribution of the y = x^ 

dx 
dy 



p{x'^) = p{x) 



Ao _i_ 
oc — y 
2 " 



(30) 



(31) 



Fig. 8: Levy fiight simulation results for the distribution of 
-^ for a = 1.5 and K^^ ^ r{l- ^), A = I, t = 10^ (left) 

and A — 100, t = 10^ (right), and theoretical predictions 
according to Eq. (l33t (solid curves). 



Hence, this example illustrates that the TAMSD of 
Levy Flights with jump length distributions without sec- 
ond moment is a one-sided Levy density lacking the first 
moment. As will be shown later, the TAMSD distribution 
for the corresponding Levy walk is competely different de- 
spite the similarity in the central part of the propagator, 
see Fig. [HI 



4 Enhanced diffusion regime 

In this section we consider the regime of sojourn times 
distributed according to a PDF with existing mean (r), 
but diverging second moment, i.e. Eq. ^ with 1 < a < 2. 
The expansion in Laplace domain is hence 



Vi(w) ~ 1- {t)u + Au°'. 



(34) 



Note that the transition from x to the positive valued 
y results in a factor 2 in the normalization. The large 



In our case Eq. (0) the average sojourn time is (t) = 

a/(a-l) and A^ \r{l ~ a)\. 
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4.1 Particle position distribution 

Again, the particle position is given by the integral over 
the velocities vq (i_(- — i- ) • The MSD is well known [3D] : 



\ ^'^ - (r) r(4-a) ' 



(35) 



as well as the propagator which is given by a two-sided 
Levy-distribution for the central part where \x\ <ti vqI, 



p{x,t) 



1 



(i^„t) 



l/a 



La,0 



(K^t) 



l/a 



(36) 



where La^ik) — exp [— fc"] and 

vlA{a - 1) 



i^. 



{T)r{A-ay 



(37) 



stationarity at large ti and obeys the limiting distribution 

m 

1 /""^ A 

ti->oo (r) J^^ (t) |7 (1 - a)\ ■* 

This first waiting time in turn has no mean and therefore 
Po dominates {viti)v{ti+ A)). Again, with Eqs. ^, (fTO)) . 
(w(ii)u(ti + Zi)) = VQpo{ti,ti + A) holds for the velocity 
correlation function for ti and A both large. The po for 
the present process is well known, so that following e.g. 
the procedure presented in [30] we have 



{v{ti)v{ti + A)) = 

h 

h + A 



v^A 



(a-l)(r)|r(l-a; 

l-Q 



rtj"" X 



h 



ti + A 



(38) 



In the equilibrated regime (or stationary state) ii ^ A, 
{v{ti)v{ti + A)) becomes independent of time ii so that 




{V{tl)v{h + A))^g 



v'oA 



1 1 — a 



Fig. 9: For the Levy walk with 1 < a < 2 the central 
part of the particle distribution (vs. the scaling variable 
z = x/t^'") is a symmetric Levy distribution Eq. (|36l) 
(solid line). Simulation (histogram) for a — 3/2, t — 10^, 
uo = 1, Ka see Eq. ([57]). sample size 10^. 



Due to the finite velocity the Levy walk propagator 
exhibits a cutoff so that p{x,t) = at |a;| > v^t, simi- 
larly to other Levy walk models |2S], [3S]. Moreover, we 
expect those particles that have never changed their di- 
rection of motion to form a delta-peak at the edge of the 
cutoff ^37i. In the following we will calculate the time av- 
eraged mean squared displacement (TAMSD) of the Levy 
Walk described above. 



4.2 Ensemble average of 5'^ 

It is important to note that also in the subballistic regime 
1 < a < 2, the velocity correlation is governed by the per- 
sistence probability po{t) to stay in one state for a time t: 
For sojourn time PDFs ■(/'(''') with existing mean the corre- 
sponding forward recurrence time PDF V'/.ti ('''/) reaches 



{T){a-l)\r{l-a) 



.(39) 



By integrating Eq. ([55)1 as in Eq. ([5|) we obtain the position 
autocorrelation 



{x{h)x{h+A)) 



Avl 



{r)nA~a) ' 



i?-" X 



{yf-^ + (1 + yf-" + (« - 3) (1 + 2/)'-" + ^[40) 



where y — A/ti. Our simulations have shown that this 
estimation reproduces the large ii behavior of the position 
correlation {x{ti)x{t2)) quite well. For Z\ — ;> the behavior 
of the MSD Eq. ^ is reproduced. 

Note that in the subballistic case the MSD for a pro- 
cess starting with the beginning of the measurement dif- 
fers from the MSD for a process that started a long time 
before the beginning of the measurement time io, as was 
found earlier in the context of a stochastic collision model 
[38] . This behavior is due to the predominant role of the 
persistence probability po in the correlation functions dis- 
cussed above. In what follows the MSD for the process 
that started long before to will be called the equilibrium 
MSD (x^) , alluding to the fact that this process has no 
memory of its starting time. The situation is sketched in 
Fig. [TU] It is important to introduce the equilibrium MSD 
at this point for the following reason: Since the time aver- 
aging procedure comprises averaging over all continuously 
shifted time lags, and not only over those starting at a 
switching event, there is an inherent averaging over dis- 
order. The equilibrium MSD accounts for this averaging 
over disorder and is therefore the natural ensemble aver- 
aged quantity to later compare the time averaged MSD to. 
Note also that such a definition of an equilibrium MSD is 
not possible in the ballistic case since there a stationary 
state does not exist - the respective MSD would never be- 
come independent of the time difference between start of 
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T1- 



^ T2— «T3-H- U _^ 



Numerical simulations also indicate that these two aver- 
ages appear to differ by a factor (Fig. [TT|) . Further nu- 
merical evidence for this behavior was found in [23], |24| . 
Especially for small a the convergence of (S^) is extremely 
slow. Moreover, in our simulations the large time behav- 
ior of the (x^) is not represented very accurately at the 
corresponding relatively small times (up to 10"*). 



fo 



^ T2— hT3.|_ U _^| 



Fig. 10: Sketch of the process starting at the beginning 
of the measurement tg (upper panel) and the equilibrium 
process starting at t — —tmi with tmi <C (r). The mea- 
surement begins at a time Iq between two renewal events 
(lower panel). Renewal events are indicated by black ticks, 
the start of the measurement to is marked by a cross. 



the process in the past and the actual start of the mea- 
surement. 

Using Eq. ^ and (x^)^^ = 2 /^ dh J* dt2{v{h)v{t2))cq 
we obtain the known equilibrium MSD |38] 



A 
V 



3- 



▲ A 



100 200 



500 



1000 
A 



2000 5000 1x10'' 



Fig. 11: Ratio (<5^)/(a;^) for a — 1.25 (triangles), a = 

1.5 (diamonds) and a = 1.7 (circles); t — 10^. The theory 
predicts 4, 2 and 1.43 for this ratio (solid lines). 



{^\ 



A 



j-3 — a 



(r)r(4-a) 



(41) 



For the TAMSD we use again the definition Eq. ([T]) and 
write Eq. ([7]) for the respective ensemble average. To ob- 
tain a description of the ensemble averaged TAMSD, we 
insert the autocorrelation function Eq. (HUl) and the MSD 
Eq. ([33]) into Eq. ©. We thus have 



2A 



^ ' {t^A){r)r{5-a) 
- (t - Af'" + t^-" - A'^-°' 
-^(4 - aYA^-"" - (4 - a)t3-"Z\] 



(42) 



which becomes in leading orders by expansion for A/t 
small: 



(<52) 



2Avl 



(T)r(4 - a) 



A 



3-Q 



Avl 



{T)r{s 



-AH^ 



(43) 



Note that this result for ((5^) complies with the time de- 
pendence of the equilibrium ensemble average {x^) , Eq. 
(PT|) . and hence differs from the MSD p5p by lacking a 
factor: 

'P 1 
lim -r-^ = r. (44) 



t^ca (x^) 



1 



4.3 Fluctuations of the TAMSD 

The TAMSDs of trajectories measured up to a certain ob- 
servation time appear to be distributed. Fig. [T^ shows the 
TAMSD evolution of some sample trajectories. The fluc- 
tuations of the TAMSDs decrease with increasing total 
measurement time. However, since in experiments the ob- 
servation time is always finite, these fluctuations may play 
a role in practice. In our simulations for example we find 
large fluctuations among the 5'^ for a = 1.25 and A = 100, 
t= 10^ (see Fig. [12]). 

Although the propagators of the Levy walk and flight 
look very similar in the central part \x\ < vot, unlike the 
flight case simulations suggest that the TAMSD distribu- 
tion for the Levy walk cannot be expressed by a Levy dis- 
tribution of stability index a/2. For the Ergodicity Break- 
ing (EB) parameter of the subballistic Levy Walk regime. 



EB 



lim 

i— >-oo 



{Pj 



-{5^ 



5-^ 



(45) 



we find a steady decay with t, which is yet very slow with 
an exponent of roughly (1 — a), to the value zero indicat- 
ing that the width of the distribution tends to zero (Figs. 
[131 for a = 1.5 , 1.25). Hence, the TAMSDs do not remain 
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Fig. 12: TAMSDs of particle trajectories in dependence on 
the lag time A, t = 10^; a = 5/4. For finite measurement 
times t fluctuations in J^ are observed. Larger times t and 
smaller lag times A result in smaller fluctuations. 



distributed in the limit of very long times in the subbal- 
listic regime. Such a very slow decay to ergodic behavior 
is not a distinct feature of the enhanced phase of the Levy 
walk model but can also be found for completely different 
ergodic systems such as relaxation of confined fractional 
Brownian motion 1391. 
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Fig. 13; EB-parameter Eq. (gS]), for a = 3/2 (left) and 5/4 
(right). Blue diamonds indicate A — 100, green triangles 
A ~ 5000, the grey solid line indicates the slope (1 — a). 



In contrast to the flight case, the width of the S'^~ 
distribution for the Levy walk at finite times always exists 
due to finite velocity. Simulations suggest that it increases 
with the lag time Z\ as Z\^. The width of the TAMSD dis- 
tribution hence appears to have the same Z\-dependence 
as a ballistic motion. However, it is not quite clear whether 
this dependence represents the large time behavior of the 
distribution of the TAMSDs, or whether it is an artefact 
of the ballistic tails of the propagator due to the extremely 
long transients. 



5 Summary 

In this article we have shown that the shifted time aver- 
aged MSD of the ballistic L'evy walk is described by the 
Mittag-Leffler distribution, similar to the distribution of 
the TAMSD in the sub-diffusive continuous time random 
walk (CTRW). This distribution describes the fluctuations 



of the time averages and is universal. The TAMSD aver- 
aged over an ensemble of trajectories {S"^) is not equal 
to the ensemble average (x^) as already pointed out by 
Akimoto [53]. Interestingly (S^) — {vqS'^) exhibits two be- 
haviors valid for A < 1 and Z\ > 1, and it would be in- 
teresting to see if a similar cross-over takes place in other 
models such as the sub-diffusive CTRW. For Levy flights 
the TAMSDs are random with a PDF given by the one 
sided Levy PDF, which is in agreement with rigorous re- 
sults (though note that our coefficients are different than 
those reported in [3S], possibly due to a typo). In the en- 
hanced diffusion regime, the PDF of the particle position 
of the Levy walk is similar to the Levy flight case, at least 
in its center. However, the TAMSD of the two models is 
vastly different, and for Levy walks no fluctuations are 
found for S'^. This indicates that the TAMSD is controlled 
by rare events, since the tails of the mentioned distribu- 
tions are where one flnds differences between the models. 
Thus taking into consideration finite velocity (like in the 
Levy walk model) is crucial for our understanding of the 
ergodic properties of these processes. In the future it might 
be worth while checking the time averages of lower order 
moments, since they might exhibit behavior very differ- 
ent the second moment considered here. We note that for 
finite times the fiuctuations of TAMSDs are large. Con- 
sequently, in the laboratory where experiments are made 
for finite time the process may seem non ergodic, but this 
is only a finite time effect. Moreover, in this sub-ballistic 
case ((5^) is equal to the equilibrium MSD {x'^)cq, but not 
to (x^). If one wishes to compare time and ensemble av- 
erages, the conclusion on equality of these two averages 
will depend on how the ensemble is prepared. To attain 
ergodicity we need to start the ensemble in a stationary 
state, which is not so surprising. The point is that for 
normal processes, e.g. the case where all moments of the 
waiting time PDF iP{t) exists, it does not matter how we 
start the process, in the long time limit the time and en- 
semble average procedures are all identical. In that sense 
the Levy walk, in the enhanced regime, is unique. A sim- 
ilar effect might be found also for sub-diffusive CTRW, 
for the case of finite average waiting time, but with an 
infinite variance, but that is left for future work. Further 
time averaged drifts, when a bias is present, also exhibit 
interesting ergodic features and Einstein relations, as we 
discussed recently in |25| . 

This work was supported by the Israel Science Foundation. 



A Correlation functions 

For the derivation of the velocity correlation function we follow 
[30| . Hence, for the time variables f i , A with t2 — ti — A and 
going to Laplace domain with respect to A we have 



Pn{tl,UA) = 






, i-i'(u^) 







n > 1 



(46) 



where i/'/,ti is the forward recurrence time PDF, i.e. the PDF 
of the time it takes to encounter the next event after a given 
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time ii wiiicii in double-Laplace domain reads 

1 iP{ua) ~ i'{ui) 



Ipf.uiiuA) 



1 - Tp{ui) Wl - IJ-A 



(47) 



and repeated integration by parts we find 

,,,,., 2 sin Tva 
[x{ti)x{t2)) = Vo X 



dsi / B [ — ; a, 1 — a ) ds2 
Js, \*2 



Here ua and mi are the Laplace variables conjugate to A and 
ti, respectively. Inserting ()46[) and H47[) back into ((9| we get 
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■"0 
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l-i>(ui) ui-UA ua(i + 4'{ua)) 



2 1 

^ Vo — 
UA 



1pf,ui{uA) —, 

^1 1 + ^P{ua) 



(48) 



where we denote the Laplace transformation by £ 
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2 sm 7ra 

Vo X 
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1 + Q, 1 — a 
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^:t?B ( -, 1 + Q, 1 - Q ) + l-tlB (1, -1 + a, 1 - a) 



t2' 



+ |(l-a)i? 



(51) 



which with 5(1, — l + a, 1 — a) =0 for < a < 1 finally yields 
Eq. JISl). 



A.l Ballistic phase 



Let us now specify the sojourn time distribution at large times, 
in Laplace domain '4'(u) — 1 — Au" , < 1 < a which yields 



Cti,A {{v{ti)v{ti + A))\ui,ua} = 

1 1 ^{uA)~i^{ux) 



A. 2 Subballistic phase 

Here the Laplace transform of the sojourn time of the particle 
in a velocity state for small u reads 'tp{u) = 1 — {t)u + Au" 
with (r) = a. I a — 1 and A = \r(\ — q)|. For t\,A large, i.e. 
u\,UA — >■ in Laplace domain, Eq. (|48|) becomes again 
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UAU\ l-?/i(ui) (u\~ Ua)ua 



2 

uo 



i>(UA) - i>(u\) 



UAUl l_^(i(^) {ui-Ua)ua 



2 

fo 



1 



UAUl 



i>f,ui{uA) 

UA 



(49) 
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■ VoVo\Ui,Ua) = Vo 



and after Laplace inversion 
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Laplace inversion yields 
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(50) 
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-i)|r(i- 
1 


-a)\ 



which we will use in the following since we know the scaling 
form of V'/,ti {^) for large ii and A due to Dynkin's theorem 
Eq. HTl, leading to (fT2)l . 

Inserting Eq. (|12|l into Eq. ^ and using the definition of 
the incomplete Beta function B{y; a, b) — L" duu°~^{l — u)*"^^ 



(r)(Q-l)|r(l-a 
Inserting this again into Eq. (|5J gives 

(a::(ti)a:(t2)> = 
A 1 



(52) 



[A'-°' - {ti + Ay-"] 

[(i2-ii)^-°-4-°]. (53) 
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= uo- 



(r)r(4-a) 
which finally yields Eq. ((40] 
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B Ensemble averaged TAMSD, ballistic phase 



Inserting Eq. (jlSp and Eq. (jldp into Eq. (0 and again using in- 
tegration by parts and the integral definition of tlie incomplete 
Beta function yields 
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For the small A expansion Eq. (|16|l we used the identity 
B{l,b,a] 



B[^-^^,a,b 



B\ ^,b,a 



function for small arguments y 



B{y,a,b) =2/°^ 



(1-^). 
n\{a + n) 



(57) 



where (c)^ = r{c + n)/r(c) is the Pochhammer symbol. 
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